library(raster)
library(rgdal)
raster_filepath <- system.file("raster/srtm.tif", package = "spDataLarge")
new_raster <- raster(raster_filepath)
new_raster
# class : RasterLayer
# dimensions : 457, 465, 212505 (nrow, ncol, ncell)
# resolution : 0.0008333333, 0.0008333333 (x, y)
# extent : -113.2396, -112.8521, 37.13208, 37.51292 (xmin, xmax, ymin, ymax)
# crs : +proj=longlat +datum=WGS84 +no_defs
# source : srtm.tif
# names : srtm
# values : 1024, 2892 (min, max)Spatial Information Analysis
GIS
Code
R
공간정보분석
Geocomputation with R
Chapter 1 : R을 이용한 공간정보 분석
1. GIS와 공간정보 데이터
1.1 공간정보와 GIS
공간정보란 ? 사람들이 생활하고 있는 공간 상에서 사건이나 사물에 대한 위치를 나타내는 정보
위치를 나타내는 정보는 (1) 위치를 표현하는 정보 (2) 해당 위치에 나타나는 특성에 대한 정보
위치를 표현하는 정보 : 공간 상에서 사건이나 사물의 위치가 어디에 있는지를 나타내는 정보
- ex) 주소, 위경도, x,y 좌표 등
해당 위치에 나타나는 특성에 대한 정보 : 특정 위치에 있는 사건이나 사물을 설명하는 정보
- ex) 학교, 회사, 학생 수 , 교사 수, 사고 건 수, 사고 유형 등
지리정보 시스템(Geographic Information System) : 공간정보데이터를 처리, 가공하여 새로운 정보를 도출하는 일련의 과정 또는 기법
- ex) 교통사고 데이터 분석 (TIMS)
공간정보를 이용하여 GIS 분석을 수행하기 위한 소프트웨어
- 전용 소프트웨어
- ArcGIS : 전문적인 공간정보의 처리와 분석 가능, 고가(유료)
- QGIS : 오픈소스 GIS 소프트웨어, 최근 많은 분야에서 GIS 소프트웨어로 활용
- 오픈소스 소프트웨어
- R 소프트웨어 : 오픈소스 기반의 통계 프로그램, 공간정보의 처리와 븐석에도 강력한 기능
- Python 소프트웨어 : 배우기 쉽고, 강력한 프로그래밍 언어, 공간정보를 다루는데 유용한 라이브러리가 개발
- 전용 소프트웨어
1.2 공간정보 데이터
- 위치정보와 속성정보로 구분
- 위치정보
- 좌표체계를 이용한 위치정보
- 지리좌표계에서 이용하는 경도와 위도로 표현 ex) 경위도좌표
- 수학적으로 X좌표와 Y좌표로 위치 정보를 표현 ex) 평면직각좌표(지도좌표)
- 공간정보 데이터의 위치정보 표현 방식
- 벡터 (점, 선, 면)
- 래스터 (일정한 격자 또는 화소)
- 좌표체계를 이용한 위치정보
- 속성정보
- 주어진 위치에 있는 사건이나 사물에 대한 자료
- 위치정보
1.3 공간정보 좌표체계
- 지리좌표체계 : 경도와 위도로 위치를 표현하는 지리좌표체계
- 투영좌표체계 : 지도투영법을 적용하여 둥근 지구를 평면으로 변환한 후, 직각좌표체계를 이용하여 x좌표와 y좌표의 직각좌표체계로 위치를 표현
- 원통도법, 원추도법, 평면도법이 있음.
- UTM 좌표체계, TM 좌표계, UTM-K 좌표계
- 우리나라는 ITRF2000 지구중심좌표계를 따르고 타원체로는 GRS80 타원체를 적용
1.4 공간정보 파일
- shapefile
- .shp : 공간정보(점, 선, 다각형)
- .shx : geometry와 속성 정보 연결
- .dbf : 속성정보
- .drj : 좌표계 정보 저장
- .sbn : 위치 정보 저장
- geojson : json 또는 xml 파일 포맷 필요요
Chapter2 : Geographic data in R
패키지
sf: 지리 공간 벡터 데이터(vector data) 분석을 위한 패키지raster: 지리 공간 레스터 데이터(raster data)를 처리 및 분석하는데 사용spData: 37개의 지리 공간 데이터셋이 내장spDataLarge: 지리공간 데이터 샘플을 내장
vignetee(package = " "): 설치된 모든 패키지에 대한 이용가능한 모든 목록을 출력st_as_sf(): st 데이터를 sf로 변환하는 함수st_centroid: 폴리곤의 중심점을 계산하는 함수plot 함수 위에 다른 지도 층을 추가 :
plot()함수 안에add = TRUE사용
2.1 Simple feature geometries (sfg)
st_point(): A pointst_linestring(): A linestringst_polygon(): A polygonst_multipoint(): A multipointst_multilinestring(): A multilinestringst_multipolygon(): A multipolygonst_geometrycollection(): A geometry collection
2.2 Simple feature columns (sfc)
st_sfc(): 두 개의 지리특성(feature)을 하나의 칼럼 객체로 합치는 함수st_geometry_type(): 기하유형을 확인st_crs(): 특정 CRS를 지정- 특정 CRS를 지정하기 위해
epsg(SRID)또는proj4string속성을 사용
- 특정 CRS를 지정하기 위해
epsg코드- 장점 : 짧아서 기억하기 쉬움
sfc객체 내의 모든 geometries는 동일한 CRS를 가져야 함.- EPSG : 4326 : GPS가 사용하는 좌표계
proj4string정의- 장점 : 투사 유형이나 datum, 타원체 등의 다른 모수들을 구체화할 수 있는 유연성이 있음
- 단점 : 사용자가 구체화를 해야하므로 길고 복잡하며 기억하기 어려움
st_sf(): sfc와 class sf의 객체들을 하나로 통합
dim(): 행, 열, 층의 수ncell(): 셀의 수res(): 해상도extent(): 경계값crs(): 좌표계inMemory(): 래스터 데이터가 메모리에 저장되어 있는지(논리값 출력)
2.3 Raster classes
- RasterLayer class
- RasterBrick Class
- RasterStack class
RasterLayer : 한 개의 층으로 구성되어 있는 래스터
RasterBrick : 여러개의 층으로 구성되어 있는 래스터
- 단일 다중 스펙트럼 위성 파일, 메모리의 단일 다층 객체의 형태
brick()함수를 사용하여 다층 래스터 파일을 로드
RasterStack : 여러개의 층으로 구성되어 있는 래스터
nlayers(): 래스터 데이터의 층의 수
2.3.1 RasterBrick과 RasterStack의 차이
- RasterBrick : 동일한 복수 개의 RasterLayer 층으로 구성
- RasterStack : 여러 개의 RasterLayer과 RasterBrick 객체가 혼합
2.3.2 언제 어떤 래스터 클래스를 사용하는 것이 좋은가 ?
- RasterBrick : 하나의 다층 래스터 파일이나 객체를 처리
- RasterStack : 여러 개의 래스터 파일들이나 여러 종류의 래스터 클래스를 한꺼번에 연걸해서 연산하고 처리
2.4 CRS(Coordinate Reference Systems)
지리 좌표계
- 위도와 경도를 이용해 지구 표면의 위치를 정의
- 미터가 아니라, 각도로 거리 측정
- 타원 표면, 구면 표면
- WGS84
투영(투사) 좌표계
- 암묵적으로 “평평한 표면” 위의 데카르트 좌표 기반 -> 왜곡 발생
- 원점, x축, y축
- 미터와 같은 선형 측정 단위
- 평면, 원뿔, 원통의 3가지 투영 유형
st_set_crs(): 좌표계가 비어있거나 잘못 입력되어 있는 경우에 좌표계를 설정st_transform(): 투영 데이터 변환st_area(): 벡터 데이터의 면적 계산 -> [m^2] 단위가 같이 반환좌표계 설정할 때,
- 벡터 데이터 :
epsg코드나proj4string정의 모두 사용 가능 - 래스터 데이터 :
proj4string정의만 사용
- 벡터 데이터 :
Chapter3 : Attribute data operations
3.1 벡터 데이터에서 속성 정보를 가져오는 방법
sf객체에서 속성 정보만 가져오기 :st_drop_geometry()
- Base R 구문으로 벡터 데이터 속성 정보의 행과 열 가져오기
- dplyr로 벡터 데이터 속성 정보의 행과 열 가져오기
- 한 개 컬럼만 가져온 결과를 벡터로 반환하기
3.1.1. sf 객체에서 속성 정보만 가져오기 : st_drop_geometry()
- 지리공간
sf객체는 항상 점, 선, 면 등의 지리기하 데이터를 리스트로 가지고 있는 geometry 칼럼이 항상 따라다님 sf객체로부터 이 geometry 칼럼을 제거하고 나머지 속성 정보만으로 Dataframe을 만들고 싶다면sf패키지의st_drop_geometry()를 사용- geometry 칼럼의 경우 지리기하 점, 선, 면 등의 리스트 정보를 가지고 있어 메모리 점유가 크기때문에, 사용할 필요가 없다면 geometry 칼럼을 제거하고 속성 정보만으로 Dataframe으로 만들어서 분석을 진행하는게 좋음
3.1.2. Base R 구문으로 벡터 데이터 속성 정보의 행과 열 가져오기
R Dataframe에서 i행과 j열을 가져올 때 :
df[i, j],subset(),$을 사용- i행과 j열 위치를 지정 ex)
world[1:6, ]
- i행과 j열 위치를 지정 ex)
- j행의 이름을 이용 ex)
world[, c("name_long", "lifeExp")]
- j행의 이름을 이용 ex)
- 논리 벡터를 사용해서 i행의 부분집합 ex)
sel_area <- world$area_km2 < 10000
- 논리 벡터를 사용해서 i행의 부분집합 ex)
3.1.3. dplyr로 벡터 데이터 속성 정보의 행과 열 가져오기
dplyr 패키지에서는 체인(%>%)으로 파이프 연산자를 사용하여 가독성이 좋고, 속도가 빠름
select()함수를 사용하여 특정 열 선택
select(sf, name)select(sf, name1:name2)select(sf, position)ex) select(world, 2, 7)select(sf, -name)select(sf, name_new = name_old): 열 선택하여 이름 변경select(sf, contain(string)): 특정 문자열을 포함한 칼럼을 선택contain(),starts_with(),ends_with(),matches(),num_range()
filter()함수를 사용하여 조건을 만족하는 특정 행 추출
subset()함수와 동일한 기능
aggregate()함수를 사용하여 지리 벡터 데이터의 속성 정보를 그룹별로 집계
aggregate(x ~ group, FUN, data, ...)- data.frame을 반환하며, 집계된 결과에 지리 기하(geometry) 정보는 없음
- world[‘pop’]은 “sf” 객체이기 때문에 집계 결과가 “sf” 객체로 반환
- world$pop은 숫자형 벡터이므로
aggregate()함수를 적용하면 집계 결과가 “data.frame”으로 반환
summarize(),group_by()함수를 이용한 지리벡터 데이터의 속성 정보를 그룹별로 집계
group_by(): 기준이 되는 그룹을 지정summarize(): 다양한 집계 함수를 사용sum(),n(): 합계와 개수 집계top_n(): 상위 n개 추출arrange(): 오름차순 정렬,desc()를 사용하면 내림차순 정렬st_drop_geometry(): geometry 열 제거
3.1.4. 두 개의 지리 벡터 데이터 테이블을 Join하기(결합)
R의 sf클래스 객체인 지리공간 벡터 데이터를 dplyr의 함수를 사용해서 두 테이블을 join하면 속성과 함께 지리공간 geometry 칼럼과 정보도 join된 후의 테이블에 자동으로 그대로 따라감
left_join()시 key variable이 있어야 함- 두 데이터 셋에 같은 이름을 가지는 변수가 없는 경우
- 하나의 key variable의 이름을 바꿔서 통일시켜줌
by를 사용하여 결합변수를 지정
- 두 데이터 셋에 같은 이름을 가지는 변수가 없는 경우
# coffee_data의 name_long변수 이름을 nm으로 변경
coffee_renamed <- rename(coffee_data, nm = name_long)
# by 사용하여 결합 변수를 지정하여 다른이름변수를 기준으로 조인하기
world_coffee1 <- left_join(world, coffee_renamed, by = c(name_long = "nm"))inner_join()함수를 사용하면 겹치는 행만 추출setdiff(): 일치하지 않는 행 추출grepl(): 텍스트 찾는 함수 (논리값으로 출력)grep(): 텍스트 찾는 함수 (행 번호 출력)
3.1.5. 지리 벡터 데이터의 새로운 속성 만들기 및 지리정보 제거하기
- dplyr로 지리 벡터 데이터에 새로운 속성 만들기
mutate(): 기존 데이터 셋에 새로 만든 변수(열) 추가transmute(): 기존의 열은 모두 제거하고 새로 만든 열과 지리기하 geometry열만을 반환
- tidyr로 지리 벡터 데이터의 기존 속성을 합치거나 분리하기
unite(data, 병합 열, sep = "_", remove = TRUE): 기존 속성 열을 합쳐서 새로운 속성 열을 만듦- remove = TRUE를 설정해주면 기존의 합치려는 두 개의 열은 제거되고, 새로 만들어진 열만 남음
separate(): 기존에 존재하는 열을 구분자를 기준으로 두 개의 열로 분리
world_unite <- world %>%
unite("con_reg", continent:region_un, sep = ":", remove = TRUE)
names(world_unite)
# "iso_a2" "name_long" "con_reg" "subregion" "type"
# "area_km2" "pop" "lifeExp" "gdpPercap" "geom"
world_separate <- world_unite %>%
separate(con_reg, c("continent", "region_un"), sep = ":")
names(world_separate)
# "iso_a2" "name_long" "continent" "region_un" "subregion" "type"
# "area_km2" "pop" "lifeExp" "gdpPercap" "geom" - dplyr로 지리 벡터 데이터의 속성 이름 바꾸기
rename(data, new_name = old_name): 특정 속성 변수 이름 변경setNames(object = nm, nm): 여러개의 속성 칼럼을 한꺼번에 변경 또는 부여
world %>% rename(name = name_long)
new_names <- c("i", "n", "c", "r", "s", "t", "a", "p", "l", "gP", "geom")
world %>% setNames(new_names)3.2 래스터 객체 조작
- 래스터 객체의 데이터 속성은 숫자형(numeric), 정수형(integer), 논리형(logical), 요인형(factor) 데이터를 지원하며, 문자형(character)은 지원하지 않음
- 문자형으로 이루어진 범주형 변수 값을 가지고 래스터 객체의 속성을 만들고 싶으면
- 문자형을 요인형으로 변환(또는 논리형으로 변환) ->
factor()함수 사용
- 문자형을 요인형으로 변환(또는 논리형으로 변환) ->
- 요인형 값을 속성 값으로 하여 래스터 객체를 만듦
- 래스터 객체의 모든 값을 추출하거나 전체 행을 추출 :
values(),getValues()
Chapter4 : Spatial data operations
4.1 벡터 데이터의 공간 연산(Spatial operations on vector data)
4.1.1 공간 부분 집합(Spatial subsetting)
st_intersects(): 공간 부분집합 추출(교집합)
4.1.2 위상 관계(Topological relations)
st_intersects(): 공간적으로 관련이 있는 객체를 출력st_disjoint(): 공간적으로 관련되지 않은 객체만 반환st_within(): 공간적으로 완전히 객체 내부에 있는 객체들만 출력st_touches(): 공간적으로 테두리에 있는 객체들만 출력st_is_within_distance(): 공간적으로 주어진 거리보다 가까운 객체들을 반환sparse = FALSE매개변수를 설정하면 논리값으로 출력
st_intersects(p, a)
#> Sparse geometry binary predicate list of length 4, where the predicate
#> was `intersects'
#> 1: 1
#> 2: 1
#> 3: (empty)
#> 4: (empty)
st_intersects(p, a, sparse = FALSE)
#> [,1]
#> [1,] TRUE
#> [2,] TRUE
#> [3,] FALSE
#> [4,] FALSE
st_disjoint(p, a, sparse = FALSE)[, 1]
#> [1] FALSE FALSE TRUE TRUE
st_within(p, a, sparse = FALSE )[, 1]
#> [1] TRUE FALSE FALSE FALSE
st_touches(p, a, sparse = FALSE)[, 1]
#> [1] FALSE TRUE FALSE FALSE
sel <- st_is_within_distance(p, a, dist = 0.9) # can only return a sparse matrix
lengths(sel) > 0
#> [1] TRUE TRUE FALSE TRUE4.1.3 공간 결합(Spatial joining)
st_join(): 공간 결합 함수
random_joined = st_join(random_points, world["name_long"]) ; random_joined
#> Simple feature collection with 10 features and 1 field
#> Geometry type: POINT
#> Dimension: XY
#> Bounding box: xmin: -158.1893 ymin: -42.91501 xmax: 165.1157 ymax: 80.5408
#> Geodetic CRS: WGS 84
#> # A tibble: 10 × 2
#> geometry name_long
#> * <POINT [°]> <chr>
#> 1 (-58.98475 -21.24278) Paraguay
#> 2 (-13.05963 25.42744) Morocco
#> 3 (-158.1893 80.5408) <NA>
#> 4 (-108.9239 27.80098) Mexico
#> 5 (-9.246895 49.9822) <NA>
#> 6 (-71.62251 20.15883) <NA>
#> 7 (38.43318 -42.91501) <NA>
#> 8 (-133.1956 6.053818) <NA>
#> 9 (165.1157 38.16862) <NA>
#> 10 (16.86581 53.86485) Polandplot()의 옵션
- 기호(plotting symbols, characters) : pch
- 기호의 크기 : cex
- 선 두께 : lwd
- 선 유형 : lty
4.1.4 비접촉 결합(Non-overlapping joins)
any(): 특정 값이 포함되어 있는지 확인할 때 유용, 여기서 TRUE가 있는지 확인 가능
any(st_touches(cycle_hire, cycle_hire_osm, sparse = FALSE))
#> [1] FALSElibrary(mapview)
library(tmap)
tmap_mode("view")
tm_basemap("Stamen.Terrain") +
tm_shape(cycle_hire) +
tm_symbols(col = "red", shape = 16, size = 0.5, alpha = .5) +
tm_shape(cycle_hire_osm) +
tm_symbols(col = "blue", shape = 16, size = 0.5, alpha = .5) +
tm_tiles("Stamen.TonerLabels")st_transform(): 투영데이터로 변환을 위한 함수st_is_within_distance(): 임계 거리보다 가까운 객체들을 반환
cycle_hire_P <- st_transform(cycle_hire, 27700)
cycle_hire_osm_P <- st_transform(cycle_hire_osm, 27700)
sel <- st_is_within_distance(cycle_hire_P, cycle_hire_osm_P, dist = 20)
summary(lengths(sel) > 0)
#> Mode FALSE TRUE
#> logical 304 438st_join()을 사용하여dist인수를 추가하여 구할 수도 있음st_join()을 사용하면 조인된 결과의 행 수가 더 크다.- 이는 cycle_hire_P의 일부 자전거 대여소가 cycle_hire_osm_P와 여러개가 겹치기 때문임
- 겹치는 점에 대한 값을 집계하고 평균을 반환하여 문제를 해결 가능
z = st_join(cycle_hire_P, cycle_hire_osm_P,
join = st_is_within_distance, dist = 20)
nrow(cycle_hire) ; nrow(z)
#> [1] 742
#> [1] 762
z = z %>%
group_by(id) %>%
summarize(capacity = mean(capacity))
nrow(z) == nrow(cycle_hire)
#> [1] TRUE4.1.5 공간 데이터 집계(Spatial data aggregation)
aggregate()와group_by() %>% summarize()를 활용하여 그룹별 통계값 계산(평균, 합 등)
# aggregate() 사용
nz_avheight <- aggregate(x = nz_height, by = nz, FUN = mean)
plot(nz_avheight[2])
# group_by() %>% summarize() 사용
nz_avheight2 <- nz %>%
st_join(nz_height) %>%
group_by(Name) %>%
summarize(elevation = mean(elevation, na.rm = TRUE))
plot(nz_avheight2[2])
st_interpolate_aw(): 면적의 크기에 비례하게 계산(면적 가중 공간 보간)
sum(incongruent$value)
#> [1] 45.41184
agg_aw = st_interpolate_aw(incongruent[, "value"],
aggregating_zones,
extensive = TRUE)
#> Warning in st_interpolate_aw.sf(incongruent[, "value"], aggregating_zones, :
#> st_interpolate_aw assumes attributes are constant or uniform over areas of x
agg_aw$value
#> [1] 19.61613 25.668724.1.6 거리 관계 (Distance relations)
- 위상 관계는 binary인 반면 거리 관계는 연속적임
st_distance(): 두 객체 사이의 거리 계산
nz_heighest <- nz_height %>% top_n(n = 1, wt = elevation)
canterbury_centroid <- st_centroid(canterbury)
#> Warning: st_centroid assumes attributes are constant over geometries
st_distance(nz_heighest, canterbury_centroid)
#> Units: [m]
#> [,1]
#> [1,] 115540
co <- filter(nz, grepl("Canter|Otag", Name))
st_distance(nz_height[1:3, ], co)
#> Units: [m]
#> [,1] [,2]
#> [1,] 123537.16 15497.72
#> [2,] 94282.77 0.00
#> [3,] 93018.56 0.00
plot(st_geometry(co)[2])
plot(st_geometry(nz_height)[2:3], add = TRUE)
4.2 래스터 데이터의 공간 연산(Spatial operations on raster data)
4.2.1 공간 부분 집합(Spatial subsetting)
cellFromXY()orraster::extract(): 좌표값을 Cell ID로 변환
id = cellFromXY(elev, xy = matrix(c(0.1, 0.1), ncol = 2))
elev[id]
#> elev
#> 1 16
terra::extract(elev, matrix(c(0.1, 0.1), ncol = 2))
#> elev
#> 1 16
clip = rast(xmin = 0.9, xmax = 1.8, ymin = -0.45, ymax = 0.45,
res = 0.3, vals = rep(1, 9))
elev[clip]
#> elev
#> 1 18
#> 2 24
terra::extract(elev, ext(clip))- operator는 raster의 다양한 inputs을 받고,
drop=FALSE로 설정했을 때, raster 객체를 반환
elev[1:2]
#> elev
#> 1 1
#> 2 2
elev[2, 1:2]
#> elev
#> 1 7
#> 2 8
elev[1:2, drop = FALSE] # spatial subsetting with cell IDs
#> class : SpatRaster
#> dimensions : 1, 2, 1 (nrow, ncol, nlyr)
#> resolution : 0.5, 0.5 (x, y)
#> extent : -1.5, -0.5, 1, 1.5 (xmin, xmax, ymin, ymax)
#> coord. ref. : lon/lat WGS 84 (EPSG:4326)
#> source(s) : memory
#> name : elev
#> min value : 1
#> max value : 2
elev[2, 1:2, drop = FALSE] # spatial subsetting by row,column indices
#> class : SpatRaster
#> dimensions : 1, 2, 1 (nrow, ncol, nlyr)
#> resolution : 0.5, 0.5 (x, y)
#> extent : -1.5, -0.5, 0.5, 1 (xmin, xmax, ymin, ymax)
#> coord. ref. : lon/lat WGS 84 (EPSG:4326)
#> source(s) : memory
#> name : elev
#> min value : 7
#> max value : 84.2.2 Local operations
elev + elev # 더하기
elev^2 # 제곱
log(elev) # 로그
elev > 5 # 논리tmap을 활용한 시각화
- tmap을 plot하기 위해서는 우선
tm_shape()로 지정해야하며,+연산자로 레이어를 추가해야함- ex)
tm_polygons(),tm_raster(),tm_borders(),tm_symbols()등
- ex)
- Interactive maps :
tmap_mode()를 사용하여"plot",과"view"모드 사용 가능 - Facet : 하나의 창에 여러 맵을 동시에 그리기
- Facet 하는 3가지 방법
- 여러변수 이름 추가
byargument oftm_facets로 공간 데이터를 나누기tmap_arrange()사용
- Facet 하는 3가지 방법
tm_basemap(): 지도를 표현할 수 있는 바탕이 되는 지도
# 1. 여러 변수 이름 추가
tmap_mode("plot")
data(World)
tm_shape(World) +
tm_polygons(c("HPI", "economy")) +
tm_facets(sync = TRUE, ncol = 2)
# 2. by argument of `tm_facets`로 공간 데이터 나누기
tmap_mode("plot")
data(NLD_muni)
NLD_muni$perc_men <- NLD_muni$pop_men / NLD_muni$population * 100
tm_shape(NLD_muni) +
tm_polygons("perc_men", palette = "RdYlBu") +
tm_facets(by = "province")
# 3. `tmap_arrange` 함수 사용 : 각각 그린다음에 배치
tmap_mode("plot")
data(NLD_muni)
tm1 <- tm_shape(NLD_muni) + tm_polygons("population", convert2density = TRUE)
tm2 <- tm_shape(NLD_muni) + tm_bubbles(size = "population")
tmap_arrange(tm1, tm2)
tmap_mode("view")
data(World, metro, rivers, land)
tm_basemap("Stamen.Watercolor") +
tm_shape(metro) + tm_bubbles(size = "pop2020", col = "red") +
tm_tiles("Stamen.TonerLabels")- Option and styles
tm_layout(): map layout 지정tm_options()내에서 설정tmap_options_diff(): default tmap options과 차이점 출력tmap_options_reset(): default tmap options으로 설정- reset을 해주지 않으면 option이 계속 설정되어있음
tmap_style(): 지도 스타일 설정
tmap_mode("plot")
tm_shape(World) +
tm_polygons("HPI") +
tm_layout(bg.color = "skyblue", inner.margins = c(0, .02, .02, .02))
tmap_options(bg.color = "black", legend.text.color = "white")
tm_shape(World) + tm_polygons("HPI", legend.title = "Happy Planet Index")
tmap_style("classic")
## tmap style set to "classic"
## other available styles are: "white", "gray", "natural", "cobalt",
## "col_blind", "albatross", "beaver", "bw", "watercolor"
tm_shape(World) +
tm_polygons("HPI", legend.title = "Happy Planet Index")
- Exporting maps
tm <- tm_shape(World) +
tm_polygons("HPI", legend.title = "Happy Planet Index")
## save an image ("plot" mode)
tmap_save(tm, filename = "./Spatial_Information_Analysis/world_map.png")
## save as stand-alone HTML file ("view" mode)
tmap_save(tm, filename = "./Spatial_Information_Analysis/world_map.html")- Quick thematic map
qtm(World, fill = "HPI", fill.pallete = "RdYlGn")
Chapter5 : Geometry operations
5.1 벡터 데이터에 대한 기하학적 연산
5.1.1 단순화(Simplification)
- 단순화는 일반적으로 더 작은 축척 지도에서 사용하기 위한 벡터 객체(선, 다각형)의 일반화를 위한 프로세스
st_simplify(): 정점을 제거하여 선을 단순화시킴dTolerance: 단위가 m이며 커질수록 더 단순화
seine_simp <- st_simplify(seine, dTolerance = 2000) # 2000m
plot(seine)
plot(seine_simp)
object.size(seine) ; object.size(seine_simp)
#> 18096 bytes 9112 bytes
- 단순화는 다각형에도 적용 가능
st_simplify()를 사용하였을 때, 영역이 겹치는 경우도 발생- rmapshaper 패키지의
ms_simplify()함수를 사용 keep_shapes = TRUE: 개체 수는 그대로 유지
us_states
us_states2163 <- st_transform(us_states, 2163)
us_states2163
us_states_simp1 <- st_simplify(us_states2163, dTolerance = 100000)
plot(us_states[1])
plot(us_states_simp1[1])
us_states2163$AREA <- as.numeric(us_states2163$AREA)
library(rmapshaper)
us_states_simp2 <- rmapshaper::ms_simplify(us_states2163, keep = 0.01,
keep_shapes = FALSE)
plot(us_states_simp2[1])
5.1.2 중심(Centroids)
- 가장 일반적으로 사용되는 중심 연산은 지리적 중심 : 공간객체의 질량 중심
st_centroid(): 지리적 중심을 생성하지만, 때때로 지리적 중심이 상위 개체의 경계를 벗어나는 경우가 발생st_point_on_surface(): 상위 개체 위에 중심이 생성
nz_centroid <- st_centroid(nz)
seine_centroid <- st_centroid(seine)
nz_pos <- st_point_on_surface(nz)
seine_pos <- st_point_on_surface(seine)
plot(st_geometry(nz), main = "nz")
plot(nz_centroid ,add=T, col="black")
plot(nz_pos ,add=T, col="red")
plot(st_geometry(seine), main = "seine")
plot(seine_centroid ,add=T, col="black")
plot(seine_pos ,add=T, col="red")
5.1.3 버퍼(Buffers)
- 버퍼 : 기하학적 특징의 주어진 거리 내 영역을 나타내는 다각형
- 지리데이터 분석에 자주 활용됨
st_buffer(): 버퍼 생성 함수, 최소 두 개의 인수가 필요함
seine_buff_5km <- st_buffer(seine, joinStyle = "ROUND", dist = 5000)
seine_buff_20km <- st_buffer(seine, dist = 20000)
plot(seine,col="black", reset = FALSE)
plot(seine_buff_5km, col=adjustcolor(1:3, alpha = 0.2), add=T)
plot(seine,col="black", reset = FALSE)
col1 <- adjustcolor("red", alpha=0.2)
col2 <- adjustcolor("blue", alpha=0.2)
col3 <- adjustcolor("green", alpha=0.2)
plot(seine_buff_20km, col=c(col1,col2,col3), add=T)
5.1.4 아핀 변환(Affine transformations)
- 왜곡되거나 잘못 투영된 지도를 기반으로 생성된 geometry를 재투영하거나 개선할 때 많은 Affine 변환이 적용
- 이동 : 맵 단위로 모든 포인트가 동일한 거리만큼 이동
nz_sfc <- st_geometry(nz)
nz_shift <- nz_sfc + c(0, 100000)
plot(nz_sfc)
plot(nz_shift,add=T, col="Red")
- 배율 조정 : 개체를 요소만큼 확대하거나 축소
- 모든 기하 도형의 토폴로지 관계를 그대로 유지하면서 원점 좌표와 관련된 모든 좌표값을 늘리거나 줄일 수 있음
- 중심점을 기준으로 기하 도형의 차이 만큼을 늘리고 0.5배 줄인 다음 다시 중심점을 더해줌
nz_centroid_sfc <- st_centroid(nz_sfc)
nz_scale <- (nz_sfc - nz_centroid_sfc) * 0.5 + nz_centroid_sfc
plot(nz_sfc)
plot(nz_scale, add=T, col="Red")
- 회전 : 2차원 좌표의 회전하기 위한 회전변환행렬
matrix(c(cos(30), sin(30), -sin(30), cos(30)), nrow = 2, ncol = 2)
#> [,1] [,2]
#> [1,] 0.1542514 0.9880316
#> [2,] -0.9880316 0.1542514
rotation <- function(a){
r = a * pi / 180 #degrees to radians
matrix(c(cos(r), sin(r), -sin(r), cos(r)), nrow = 2, ncol = 2)
}
nz_rotate <- (nz_sfc - nz_centroid_sfc) * rotation(30) + nz_centroid_sfc
plot(nz_sfc)
plot(nz_rotate, add=T, col="red")
5.1.5 클리핑(Clipping)
- 공간 클리핑은 영향을 받는 일부 형상의 지오메트리 열의 변경을 수반하는 공간 부분 집합의 한 형태
b <- st_sfc(st_point(c(0, 1)), st_point(c(1, 1))) # create 2 points
b <- st_buffer(b, dist = 1) # convert points to circles
plot(b, border = "grey")
text(x = c(-0.5, 1.5), y = 1, labels = c("x", "y"), cex = 3) # add text
st_intersection(): X∩Y (x와 y의 교집합)st_difference(): X-Y (x와 y의 차집합)st_union(): X∪Y (x와 y의 합집합)st_sym_difference(): (X∩Y)^c (드모르간의 법칙)
par(mfrow = c(2,2))
x <- b[1] ; y <- b[2]
# X ∩ Y
x_and_y <- st_intersection(x, y)
plot(b, border = "grey", main = "X ∩ Y")
plot(x_and_y, col = "lightgrey", border = "grey", add = TRUE)
text(x = c(-0.5, 1.5), y = 1, labels = c("x", "y"), cex = 3)
# X - Y
x_dif_y <- st_difference(x,y)
plot(b, border = "grey", main = "X - Y")
plot(x_dif_y, col = "lightgrey", border = "grey", add = TRUE)
text(x = c(-0.5, 1.5), y = 1, labels = c("x", "y"), cex = 3)
# X U Y
x_union_y <- st_union(x,y)
plot(b, border = "grey", main = "X U Y")
plot(x_union_y, col = "lightgrey", border = "grey", add = TRUE)
text(x = c(-0.5, 1.5), y = 1, labels = c("x", "y"), cex = 3)
# (X ∩ Y)^c
x_sdif_y <- st_sym_difference(x,y)
plot(b, border = "grey", main = "(X ∩ Y)^c")
plot(x_sdif_y, col = "lightgrey", border = "grey", add = TRUE)
text(x = c(-0.5, 1.5), y = 1, labels = c("x", "y"), cex = 3)
5.1.6 부분집합과 클리핑(Subsetting and Clipping)
- 클리핑 오브젝트는 지오메트리를 변경할 수 있지만 오브젝트의 부분 집합을 지정할 수도 있으며 클리핑/하위 설정 오브젝트와 교차하는 피쳐만 반환할 수도 있음
st_sample(): x와 y의 범위 내에서 점들의 간단한 무작위 분포를 생성
bb = st_bbox(st_union(x, y))
box = st_as_sfc(bb)
set.seed(2017)
p = st_sample(x = box, size = 10)
x_and_y = st_intersection(x, y)
plot(b, border = "grey")
plot(p, add=T)
- X와 Y 둘 다와 교차하는 점만을 반환하는 방법
## 1번째방법
p_xy1 <- p[x_and_y]
plot(p_xy1, add=T, col="red")
## 2번째방법
p_xy2 <- st_intersection(p, x_and_y)
plot(p_xy2, add=T, col="blue")
## 3번째방법
sel_p_xy <- st_intersects(p, x, sparse = FALSE)[, 1] &
st_intersects(p, y, sparse = FALSE)[, 1]
p_xy3 <- p[sel_p_xy]
plot(p_xy3, add=T, col="green")
5.1.7 공간 결합(Geometry unions)
- 미국의 49개 주의 정보를 4개 지역으로 재구분
plot(us_states[6])
## 1. aggregate함수
regions <- aggregate(x = us_states[, "total_pop_15"], by = list(us_states$REGION),
FUN = sum, na.rm = TRUE)
plot(regions[2])
## 2. group_by, summarize함수
regions2 <- us_states %>% group_by(REGION) %>%
summarize(pop = sum(total_pop_15, na.rm = TRUE))
plot(regions2[2])
- 위에서
aggregate()와summarize()가 모두 지오메트리를 결합하고st_union()을 사용하면 지오메트리만을 분해
us_west <- us_states[us_states$REGION == "West", ]
plot(us_west[6])
us_west_union <- st_union(us_west)
plot(us_west_union)
texas <- us_states[us_states$NAME == "Texas", ]
texas_union <- st_union(us_west_union, texas)
plot(texas_union)
5.1.8 유형 변환(Type transformations)
st_cast(): 지오메트리 유형을 변환
multipoint <- st_multipoint(matrix(c(1, 3, 5, 1, 3, 1), ncol = 2))
linestring <- st_cast(multipoint, "LINESTRING")
polyg <- st_cast(multipoint, "POLYGON")
plot(multipoint)
plot(linestring)
plot(polyg)
st_length(linestring) # 길이 계산
# [1] 5.656854
st_area(polyg) # 면적 계산
# [1] 4
- multilinestring : 여러 개의 linestring을 하나의 묶음으로 처리
- multilinestring은 각 선 세그먼트에 이름을 추가하거나 단일 선 길이를 계산할 수 없는 등 수행할 수 있는 작업 수가 제한됨
st_cast()함수를 사용하여 하나의 multilinestring을 세 개의 linestring로 분리
linestring_sf2 = st_cast(multilinestring_sf, "LINESTRING")
linestring_sf2
#> Simple feature collection with 3 features and 0 fields
#> Geometry type: LINESTRING
#> Dimension: XY
#> Bounding box: xmin: 1 ymin: 1 xmax: 4 ymax: 5
#> CRS: NA
#> geom
#> 1 LINESTRING (1 5, 4 3)
#> 2 LINESTRING (4 4, 4 1)
#> 3 LINESTRING (2 2, 4 2)- name과 length 추가
linestring_sf2$name <- c("Riddle Rd", "Marshall Ave", "Foulke St")
linestring_sf2$length <- st_length(linestring_sf2)
linestring_sf2
#> Simple feature collection with 3 features and 2 fields
#> Geometry type: LINESTRING
#> Dimension: XY
#> Bounding box: xmin: 1 ymin: 1 xmax: 4 ymax: 5
#> CRS: NA
#> geom name length
#> 1 LINESTRING (1 5, 4 3) Riddle Rd 3.605551
#> 2 LINESTRING (4 4, 4 1) Marshall Ave 3.000000
#> 3 LINESTRING (2 2, 4 2) Foulke St 2.000000
plot(linestring_sf2[2])
5.2 래스터 데이터에 대한 기하학적 연산
5.2.1 공간 교집합(Geometric intersections)
- 다른 공간 객체에 의해 중첩된 래스터에서 값을 추출하는 방법
- 공간 출력을 검색하기 위해 거의 동일한 부분 집합 구문(많이 겹치는 부분)을 사용
drop = FALSE를 설정하여 행렬 구조를 유지- cell 중간점이 clip과 겹치는 셀을 포함하는 래스터 개체를 반환
elev <- rast(system.file("raster/elev.tif", package = "spData"))
clip <- rast(xmin = 0.9, xmax = 1.8, ymin = -0.45, ymax = 0.45,
resolution = 0.3, vals = rep(1, 9))
plot(elev)
plot(clip, add=T)
elve_clip <- elev[clip, drop = FALSE]
plot(elve_clip)
elev_raster <- rast(system.file("raster/elev.tif", package = "spData"))
rcc <- vect(xyFromCell(elev_raster, cell = 1:ncell(elev_raster))) # 셀의 중앙점 표시
xyFromCell(elev_raster,1) # 1번 셀의 중앙점 좌표
#> x y
#> [1,] -1.25 1.25
plot(elev)
plot(rcc,add=T)
plot(clip, add=T)
5.2.2 확장과 원점(Extent and origin)
- 다른 투사 및 해상도를 가진 두 이미지를 병합하려할 때 사용
extend(): 래스터 범위 확장- 새로 추가된 행과 열은 값 매개변수의 기본값(예 : NA)를 가짐
origin(): 래스터의 원점 좌표를 반환- 래스터의 원점은 좌표(0,0)에 가장 가까운 셀 모서리
elev <- rast(system.file("raster/elev.tif", package = "spData"))
elev_2 <- extend(elev, c(1,2), snap="near") # 아래/위 1행, 좌/우 2열 확장
plot(elev)
plot(elev_2, colNA="gray")
elev_3 <- elev + elev_2
#> Error: [+] extents do not match
elev_4 <- extend(elev, elev_2)
plot(elev_4, colNA="gray")
origin(elev_4)
#> [1] 0 0
origin(elev_4) <- c(0.25, 0.25)
plot(elev_4, colNA="black", add=T)
5.2.3 감소와 증가(Aggregation and disaggregation)
- 래스터 데이터 셋은 해상도가 서로 다를 수 있음
- 해상도를 match 시키기 위해 하나의 래스터 해상도를 감소(
aggregate())시키거나 증가(disagg()) 시켜야 함
# devtools::install_github("geocompr/geocompkg")
dem <- rast(system.file("raster/dem.tif", package = "spDataLarge"))
dem_agg <- aggregate(dem, fact = 5, fun = mean)
dem_disagg <- disagg(dem_agg, fact = 5, method = "bilinear")
plot(dem)
plot(dem_agg)
plot(dem_disagg)
identical(dem, dem_disagg)
#> [1] FALSE- 새롭게 만들어지는 cell의 값을 만드는 두가지 방법
- Default method(method = “near”) : 입력 셀의 값을 모든 출력 셀에 제공
- bilinear method : 입력 이미지의 가장 가까운 4개의 픽셀 중심을 사용하여 거리에 의해 가중된 평균을 계산
5.2.4 리샘플링(Resampling)
- Resampling : 원래 그리드에서 다른 그리드로 래스터 값을 전송하는 프로세스
- 이 프로세스는 원래 래스터의 값을 가지고, 사용자 지정 해상도와 원점을 가지고 대상 래스터의 새 값을 다시 계산함
- 해상도/원점이 다른 래스터의 값을 재계산(추정)하는 방법
- Nearest neighbor : 원래 래스터의 가장 가까운 셀 값을 대상 래스터의 셀에 할당. 속도가 빠르고 일반적으로 범주형 래스터에 적합
- Bilinear interpolation(이중선형보간) : 원래 래스터에서 가장 가까운 4개의 셀의 가중 평균을 대상 1개의 셀에 할당. 연속 래스터를 위한 가장 빠른 방법
- Cubic interpolation(큐빅 보간) : 본 래스터의 가장 가까운 16개 셀의 값을 사용하여 출력 셀 값을 결정하고 3차 다항식 함수를 적용. 연속 래스터에 사용. 2선형 보간보다 더 매끄러운 표면을 만들지만, 계산적으로 까다로움
- Cubic spline interpolation(큐빅 스플라인 보간) : 원래 래스터의 가장 가까운 16개의 셀의 값을 사용하여 출력 셀 값을 결정하지만 큐빅 스플라인(3차 다항식 함수)을 적용
- Lanczos windowed sinc resampling(Lanczos 윈도우 재샘플링) : 원래 래스터의 가장 가까운 셀 36개의 값을 사용하여 출력 셀 값을 결정
summin, q1, med, q3, max, average, mode, rms
- Nearest neighbor은 범주형 래스터에 적합한 반면, 모든 방법은 연속형 래스터에 사용
resample(x, y, method = "bilinear", filename = "", ...): 리샘플링 함수
library(terra)
target_rast <- rast(xmin = 794600, xmax = 798200,
ymin = 8931800, ymax = 8935400,
resolution = 150, crs = "EPSG:32717")
target_rast
#> class : SpatRaster
#> dimensions : 24, 24, 1 (nrow, ncol, nlyr)
#> resolution : 150, 150 (x, y)
#> extent : 794600, 798200, 8931800, 8935400 (xmin, xmax, ymin, ymax)
#> coord. ref. : WGS 84 / UTM zone 17S (EPSG:32717)
plot(dem)
plot(target_rast)
"near": 셀에 가장 가까운 픽셀에서 값을 가져옴
dem_resampl_1 <- resample(dem, y = target_rast, method = "near")
plot(dem_resampl_1)
"bilinear": 네 개의 가장 가까운 셀의 가중 평균
dem_resampl_2 <- resample(dem, y = target_rast, method = "bilinear")
plot(dem_resampl_2)
"average": 각각의 새로운 셀이 중복되는 모든 입력 셀의 가중 평균
dem_resampl_3 <- resample(dem, y = target_rast, method = "average")
plot(dem_resampl_3)
Chapter 6 : Raster-vector interactions
6.1 Raster cropping(잘라내기)
- 입력 래스터 데이터 세트의 범위가 관심 영역보다 클 경우 래스터 자르기(Cropping) 및 마스킹(Masking)은 입력 데이터의 공간 범위를 통합하는 데 유용함
- 두 작업 모두 후속 분석 단계에 대한 객체 메모리 사용 및 관련 계산 리소스를 줄이고 래스터 데이터를 포함하는 매력적인 맵을 만들기 전에 필요한 전처리 단계임
- 대상 개체와 자르기 개체는 모두 동일한 투영을 가져야 함
crop(): 두 번째 인수에 대한 래스터를 잘라냄mask(): 두 번째 인수에 전달된 개체의 경계를 벗어나는 값을 NA로 설정- 대부분의 경우
crop()과mask()를 함께 사용
- 대부분의 경우
srtm <- rast(system.file("raster/srtm.tif", package = "spDataLarge"))
zion <- read_sf(system.file("vector/zion.gpkg", package = "spDataLarge"))
zion <- st_transform(zion, crs(srtm)) # zion을 srtm 좌표계랑 동일하게
plot(srtm)
plot(vect(zion),add=T)
srtm_cropped <- crop(srtm, vect(zion))
plot(srtm_cropped)
srtm_masked <- mask(srtm, vect(zion))
plot(srtm_masked)
srtm_cropped <- crop(srtm, vect(zion))
srtm_final <- mask(srtm_cropped, vect(zion))
plot(srtm_final)
updatevalue = 0: 외부의 모든 픽셀이 0으로 설정inverse = TRUE: 경계 내에 있는 것들이 마스킹
srtm_update0 <- mask(srtm, vect(zion), updatevalue = 0)
plot(srtm_update0)
srtm_inv_masked <- mask(srtm, vect(zion), inverse = TRUE)
plot(srtm_inv_masked)
tmap을 활용한 시각화
## Original / Crop / Mask / Inverse Map
library(tmap)
library(rcartocolor)
terrain_colors = carto_pal(7, "Geyser")
pz1 = tm_shape(srtm) +
tm_raster(palette = terrain_colors, legend.show = FALSE, style = "cont") +
tm_shape(zion) +
tm_borders(lwd = 2) +
tm_layout(main.title = "A. Original", inner.margins = 0)
pz2 = tm_shape(srtm_cropped) +
tm_raster(palette = terrain_colors, legend.show = FALSE, style = "cont") +
tm_shape(zion) +
tm_borders(lwd = 2) +
tm_layout(main.title = "B. Crop", inner.margins = 0)
pz3 = tm_shape(srtm_masked) +
tm_raster(palette = terrain_colors, legend.show = FALSE, style = "cont") +
tm_shape(zion) +
tm_borders(lwd = 2) +
tm_layout(main.title = "C. Mask", inner.margins = 0)
pz4 = tm_shape(srtm_inv_masked) +
tm_raster(palette = terrain_colors, legend.show = FALSE, style = "cont") +
tm_shape(zion) +
tm_borders(lwd = 2) +
tm_layout(main.title = "D. Inverse mask", inner.margins = 0)
tmap_arrange(pz1, pz2, pz3, pz4, ncol = 4, asp = NA)
6.2 Raster extraction(래스터 추출)
- 특정 위치에 있는 대상 래스터와 관련된 값을 식별하여 반환
data("zion_points", package = "spDataLarge")
elevation <-terra::extract(srtm, vect(zion_points))
zion_points <- cbind(zion_points, elevation)
plot(srtm)
plot(vect(zion),add=T)
plot(zion_points,col="black", pch = 19, cex = 0.5, add=T)
#> Warning in plot.sf(zion_points, col = "black", pch = 19, cex = 0.5, add = T):
#> ignoring all but the first attribute
st_segmentize(): 제공된 density로 line을 따라 point를 추가dfMaxLength: 최대 점의 개수
st_cast(): 추가된 point를 “POINT” 형식으로 변환
zion_transect <- cbind(c(-113.2, -112.9), c(37.45, 37.2)) %>%
st_linestring() %>%
st_sfc(crs = crs(srtm)) %>%
st_sf()
zion_transect$id <- 1:nrow(zion_transect)
zion_transect <- st_segmentize(zion_transect, dfMaxLength = 250)
zion_transect <- st_cast(zion_transect, "POINT")
#> Warning in st_cast.sf(zion_transect, "POINT"): repeating attributes for all
#> sub-geometries for which they may not be constantzion_transect <- zion_transect %>%
group_by(id) %>%
mutate(dist = st_distance(geometry)[, 1])
zion_elev <- terra::extract(srtm, vect(zion_transect))
zion_transect <- cbind(zion_transect, zion_elev)- 많은 Point들 간의 거리를 산출 : 첫번째 점들과 이후의 각각의 점들 사이의 거리 계산하기
- 횡단면의 각 점에 대한 고도값을 추출하고 이 정보를 주요 객체와 결합
tmap을 활용한 시각화
library(tmap)
library(grid)
library(ggplot2)
zion_transect_line <- cbind(c(-113.2, -112.9), c(37.45, 37.2)) %>%
st_linestring() %>%
st_sfc(crs = crs(srtm)) %>%
st_sf()
zion_transect_points <- st_cast(zion_transect, "POINT")[c(1, nrow(zion_transect)), ]
zion_transect_points$name <- c("start", "end")
rast_poly_line <- tm_shape(srtm) +
tm_raster(palette = terrain_colors, title = "Elevation (m)",
legend.show = TRUE, style = "cont") +
tm_shape(zion) +
tm_borders(lwd = 2) +
tm_shape(zion_transect_line) +
tm_lines(col = "black", lwd = 4) +
tm_shape(zion_transect_points) +
tm_text("name", bg.color = "white", bg.alpha = 0.75, auto.placement = TRUE) +
tm_layout(legend.frame = TRUE, legend.position = c("right", "top"))
rast_poly_line
plot_transect <- ggplot(zion_transect, aes(as.numeric(dist), srtm)) +
geom_line() +
labs(x = "Distance (m)", y = "Elevation (m a.s.l.)") +
theme_bw() +
# facet_wrap(~id) +
theme(plot.margin = unit(c(5.5, 15.5, 5.5, 5.5), "pt"))
plot_transect
## grid 그리기
grid.newpage() #This function erases the current device or moves to a new page.
pushViewport(viewport(layout = grid.layout(2, 2, heights = unit(c(0.25, 5), "null"))))
grid.text("A. Line extraction", vp = viewport(layout.pos.row = 1, layout.pos.col = 1))
grid.text("B. Elevation along the line", vp = viewport(layout.pos.row = 1, layout.pos.col = 2))
print(rast_poly_line, vp = viewport(layout.pos.row = 2, layout.pos.col = 1))
print(plot_transect, vp = viewport(layout.pos.row = 2, layout.pos.col = 2))
zion_srtm_values <- terra::extract(x = srtm, y = vect(zion))
group_by(zion_srtm_values, ID) %>%
summarize(across(srtm, list(min = min, mean = mean, max = max)))
#> # A tibble: 1 × 4
#> ID srtm_min srtm_mean srtm_max
#> <dbl> <int> <dbl> <int>
#> 1 1 1122 1818. 2661- 단일 영역을 특성화하거나 여러 영역을 비교하기 위해 폴리곤 당 래스터 값에 대한 요약 통계 생성
6.3 Rasterization(래스터화)
- 벡터 객체를 래스터 객체의 표현으로 변환
cycle_hire_osm <- spData::cycle_hire_osm
cycle_hire_osm_projected <- st_transform(cycle_hire_osm, "EPSG:27700")
raster_template <- rast(ext(cycle_hire_osm_projected), resolution = 1000,
crs = st_crs(cycle_hire_osm_projected)$wkt) # ext : 경계값
ch_raster1 <- rasterize(vect(cycle_hire_osm_projected), raster_template,
field = 1)
ch_raster2 <- rasterize(vect(cycle_hire_osm_projected), raster_template,
fun = "length")
ch_raster3 <- rasterize(vect(cycle_hire_osm_projected), raster_template,
field = "capacity", fun = sum)
- 폴리곤 객체를 여러 줄 문자열로 casting한 후 0.5도의 해상도로 탬플릿 래스터 생성
touches = TRUE: 경계에 해당되는 래스터만 색칠(FALSE이면 경계 내부까지)
california <- dplyr::filter(us_states, NAME == "California")
california_borders <- st_cast(california, "MULTILINESTRING")
raster_template2 <- rast(ext(california),
resolution = 0.5,
crs = st_crs(california)$wkt)
california_raster1 <-
rasterize(vect(california_borders), raster_template2,
touches = TRUE) # touches = TRUE : 경계값만
california_raster2 <-
rasterize(vect(california), raster_template2)
# with `touches = FALSE` by default, which selects only cell
6.4 Spatial Vectorization(공간 벡터화)
- 공간적으로 연속적인 래스터 데이터를 점, 선 또는 다각형과 같은 공간적으로 분리된 벡터 데이터로 변환
- 벡터화의 가장 간단한 형태는 래스터 셀의 중심부를 점으로 변환하는 것
as.points(): 모든 raster grid 셀에 대해 중심점으로 반환
elev <- rast(system.file("raster/elev.tif", package = "spData"))
elev_point <- as.points(elev) %>%
st_as_sf()
plot(elev)
plot(elev_point)
contour(): 선에 해당하는 수치 표현- 등고선의 생성 : 공간 벡터화의 또 다른 일반적인 유형은 연속적인 높이 또는 온도(등온선)의 선을 나타내는 등고선 생성
dem = rast(system.file("raster/dem.tif", package = "spDataLarge"))
cl = as.contour(dem)
plot(dem, axes = FALSE)
plot(cl, add = TRUE)
plot(dem, axes = FALSE)
contour(dem, add = T) # 수치까지 표현
as.polygons(): 래스터를 다각형으로 변환하는 것
grain <- rast(system.file("raster/grain.tif", package = "spData"))
grain_poly <- as.polygons(grain) %>%
st_as_sf()
plot(grain)
plot(grain_poly)
Chapter 7 : Reprojecting geographic data
7.1 Coordinate Reference Systems(CRS)
- CRS를 설명할 수 있는 여러가지 방법
- 단순하지만 “lon/lat 좌표”와 같이 모호할 수 있는 문장
- 공식화되었지만 지금은 구식인 proj4 strings
proj=lonlat +ellps=WGS84 +datum=WGS84 +no_defs
EPSG:4326과 같이 식별되는authority:code텍스트 문자열
- -> 3번째 방법이 가장 정확(짧고 기억하기 쉬우며 온라인에서 찾기 쉬움)
st_crs("EPSG:4326")7.2 Querying and Setting coordinate systems
- 벡터 지리 데이터 객체에서 CRS를 가져오고 설정
vector_filepath <- system.file("shapes/world.gpkg", package = "spData")
new_vector <- read_sf(vector_filepath)
st_crs(new_vector)
#> Coordinate Reference System:
#> User input: WGS 84
#> wkt:
#> GEOGCRS["WGS 84",
#> ENSEMBLE["World Geodetic System 1984 ensemble",
#> MEMBER["World Geodetic System 1984 (Transit)"],
#> MEMBER["World Geodetic System 1984 (G730)"],
#> MEMBER["World Geodetic System 1984 (G873)"],
#> MEMBER["World Geodetic System 1984 (G1150)"],
#> MEMBER["World Geodetic System 1984 (G1674)"],
#> MEMBER["World Geodetic System 1984 (G1762)"],
#> MEMBER["World Geodetic System 1984 (G2139)"],
#> ELLIPSOID["WGS 84",6378137,298.257223563,
#> LENGTHUNIT["metre",1]],
#> ENSEMBLEACCURACY[2.0]],
#> PRIMEM["Greenwich",0,
#> ANGLEUNIT["degree",0.0174532925199433]],
#> CS[ellipsoidal,2],
#> AXIS["geodetic latitude (Lat)",north,
#> ORDER[1],
#> ANGLEUNIT["degree",0.0174532925199433]],
#> AXIS["geodetic longitude (Lon)",east,
#> ORDER[2],
#> ANGLEUNIT["degree",0.0174532925199433]],
#> USAGE[
#> SCOPE["Horizontal component of 3D system."],
#> AREA["World."],
#> BBOX[-90,-180,90,180]],
#> ID["EPSG",4326]]User input: CRS식별자 (WGS 84, 입력 파일에서 가져온 EPSG:4326의 동의어)wkt: CRS에 대한 모든 관련 정보와 함께 전체 WKT 문자열을 포함input요소는 유연함(AUTHORITY:CODE(ex. EPSG:4326), CRS 이름(ex. WGS84),proj4string정의)- wkt 요소는 객체를 파일에 저장하거나 좌표 연산을 수행할 때 사용되는 WKT 표현을 저장
new_vector객체가 WGS84 타원체를 가지며, 그리니치 프라임 자오선을 사용하고, 위도와 경도의 축 순서를 사용하는 것을 볼 수 있음- 이 경우 이 CRS 사용에 적합한 영역을 설명하는 USAGE와 CRS 식별자 EPSG:4326을 가리키는 ID와 같은 추가 요소도 있음
st_crs(new_vector)$IsGeographic
#> [1] TRUE
st_crs(new_vector)$units_gdal
#> [1] "degree"
st_crs(new_vector)$srid
#> [1] "EPSG:4326"
st_crs(new_vector)$proj4string
#> [1] "+proj=longlat +datum=WGS84 +no_defs"- st_crs 함수에는 유용한 기능이 하나 있는데, 사용된 CRS에 대한 추가 정보를 검색할 수 있음.
st_crs(new_vector)$IsGeographic: CRS가 지리적 상태인지 확인st_crs(new_vector)$units_gdal: CRS 단위st_crs(new_vector)$srid: 해당 ‘SRID’ 식별자를 추출(사용 가능한 경우)st_crs(new_vector)$proj4string: proj4string 표현을 추출
st_set_crs(): CRS가 없거나 잘못 설정되어 있는 경우 CRS 설정
new_vector <- st_set_crs(new_vector, "EPSG:4326") # set CRSterra::crs(): 래스터 객체에 대한 CRS를 설정- 하지만,
crs()함수를 사용하면 좌표계는 바뀌지만 값이 바뀌지는 않음.
raster_filepath <- system.file("raster/srtm.tif", package = "spDataLarge")
my_rast <- rast(raster_filepath)
crs(my_rast)
#> [1] "GEOGCRS[\"WGS 84\",\n ENSEMBLE[\"World Geodetic System 1984 ensemble\",\n MEMBER[\"World Geodetic System 1984 (Transit)\"],\n MEMBER[\"World Geodetic System 1984 (G730)\"],\n MEMBER[\"World Geodetic System 1984 (G873)\"],\n MEMBER[\"World Geodetic System 1984 (G1150)\"],\n MEMBER[\"World Geodetic System 1984 (G1674)\"],\n MEMBER[\"World Geodetic System 1984 (G1762)\"],\n MEMBER[\"World Geodetic System 1984 (G2139)\"],\n ELLIPSOID[\"WGS 84\",6378137,298.257223563,\n LENGTHUNIT[\"metre\",1]],\n ENSEMBLEACCURACY[2.0]],\n PRIMEM[\"Greenwich\",0,\n ANGLEUNIT[\"degree\",0.0174532925199433]],\n CS[ellipsoidal,2],\n AXIS[\"geodetic latitude (Lat)\",north,\n ORDER[1],\n ANGLEUNIT[\"degree\",0.0174532925199433]],\n AXIS[\"geodetic longitude (Lon)\",east,\n ORDER[2],\n ANGLEUNIT[\"degree\",0.0174532925199433]],\n USAGE[\n SCOPE[\"Horizontal component of 3D system.\"],\n AREA[\"World.\"],\n BBOX[-90,-180,90,180]],\n ID[\"EPSG\",4326]]"
cat(crs(my_rast)) # get CRS
#> GEOGCRS["WGS 84",
#> ENSEMBLE["World Geodetic System 1984 ensemble",
#> MEMBER["World Geodetic System 1984 (Transit)"],
#> MEMBER["World Geodetic System 1984 (G730)"],
#> MEMBER["World Geodetic System 1984 (G873)"],
#> MEMBER["World Geodetic System 1984 (G1150)"],
#> MEMBER["World Geodetic System 1984 (G1674)"],
#> MEMBER["World Geodetic System 1984 (G1762)"],
#> MEMBER["World Geodetic System 1984 (G2139)"],
#> ELLIPSOID["WGS 84",6378137,298.257223563,
#> LENGTHUNIT["metre",1]],
#> ENSEMBLEACCURACY[2.0]],
#> PRIMEM["Greenwich",0,
#> ANGLEUNIT["degree",0.0174532925199433]],
#> CS[ellipsoidal,2],
#> AXIS["geodetic latitude (Lat)",north,
#> ORDER[1],
#> ANGLEUNIT["degree",0.0174532925199433]],
#> AXIS["geodetic longitude (Lon)",east,
#> ORDER[2],
#> ANGLEUNIT["degree",0.0174532925199433]],
#> USAGE[
#> SCOPE["Horizontal component of 3D system."],
#> AREA["World."],
#> BBOX[-90,-180,90,180]],
#> ID["EPSG",4326]]
crs(my_rast) <- "EPSG:26912" # set CRS
london <- data.frame(lon = -0.1, lat = 51.5) %>%
st_as_sf(coords = c("lon", "lat"))
st_is_longlat(london)
#> [1] NA
london_geo <- st_set_crs(london, "EPSG:4326")
st_is_longlat(london_geo)
#> [1] TRUE7.3 Geometry operations on projected and unprojected data
- sf는 지리 벡터 데이터에 대한 클래스와 지리 계산을 위한 중요한 하위 수준 라이브러리에 대한 일관된 명령줄 인터페이스 제공
- 구면 geometry 연산을
sf:sf_use_sf(FALSE)명령으로 끄면 버퍼는 미터와 같은 적절한 거리 단위를 대체하지 못하는 위도와 경도의 단위를 사용하기 때문에 쓸모없는 출력이 됨. - 공간 및 기하학적 연산을 수행하는 것은 경우에 따라 거의 또는 전혀 차이가 없음. (ex: 공간 부분 집합) 그러나 버퍼링과 같은 거리가 포함된 연산의 경우 (구면 지오메트리 엔진을 사용하지 않고) 좋은 결과를 보장하는 유일한 방법은 데이터의 투영된 복사본을 만들고 그에 대한 연산을 실행하는 것임.
- 그 결과 런던과 동일하지만 미터 단위의 EPSG 코드를 가진 적절한 CRS(영국 국가 그리드)에 재투사된 새로운 물체가 되었음.
- CRS의 단위가 (도가 아닌) 미터라는 사실은 이것이 투영된 CRS임을 알려줌
- 구면 geometry 연산을
london_buff_no_crs <-
st_buffer(london, dist = 1) # incorrect: no CRS
london_buff_no_crs
#> Simple feature collection with 1 feature and 0 fields
#> Geometry type: POLYGON
#> Dimension: XY
#> Bounding box: xmin: -1.1 ymin: 50.5 xmax: 0.9 ymax: 52.5
#> CRS: NA
#> geometry
#> 1 POLYGON ((0.9 51.5, 0.89862...
london_buff_s2 <-
st_buffer(london_geo, dist = 1e5) # silent use of s2 (1e5 : 10^5m = 100,000m)
london_buff_s2
#> Simple feature collection with 1 feature and 0 fields
#> Geometry type: POLYGON
#> Dimension: XY
#> Bounding box: xmin: -1.552818 ymin: 50.59609 xmax: 1.356603 ymax: 52.40393
#> Geodetic CRS: WGS 84
#> geometry
#> 1 POLYGON ((-0.3523255 52.392...
london_buff_s2_100_cells <-
st_buffer(london_geo, dist = 1e5, max_cells = 100)
london_buff_s2_100_cells
#> Simple feature collection with 1 feature and 0 fields
#> Geometry type: POLYGON
#> Dimension: XY
#> Bounding box: xmin: -1.718303 ymin: 50.51128 xmax: 1.524546 ymax: 52.53186
#> Geodetic CRS: WGS 84
#> geometry
#> 1 POLYGON ((-0.3908656 52.531...
sf::sf_use_s2(FALSE)
#> Spherical geometry (s2) switched off
london_buff_lonlat <-
st_buffer(london_geo, dist = 1) # incorrect result
#> Warning in st_buffer.sfc(st_geometry(x), dist, nQuadSegs, endCapStyle =
#> endCapStyle, : st_buffer does not correctly buffer longitude/latitude data
#> dist is assumed to be in decimal degrees (arc_degrees).
sf::sf_use_s2(TRUE)
#> Spherical geometry (s2) switched on
london_proj <- data.frame(x = 530000, y = 180000) %>%
st_as_sf(coords = 1:2, crs = "EPSG:27700")
st_crs(london_proj)
#> Coordinate Reference System:
#> User input: EPSG:27700
#> wkt:
#> PROJCRS["OSGB36 / British National Grid",
#> BASEGEOGCRS["OSGB36",
#> DATUM["Ordnance Survey of Great Britain 1936",
#> ELLIPSOID["Airy 1830",6377563.396,299.3249646,
#> LENGTHUNIT["metre",1]]],
#> PRIMEM["Greenwich",0,
#> ANGLEUNIT["degree",0.0174532925199433]],
#> ID["EPSG",4277]],
#> CONVERSION["British National Grid",
#> METHOD["Transverse Mercator",
#> ID["EPSG",9807]],
#> PARAMETER["Latitude of natural origin",49,
#> ANGLEUNIT["degree",0.0174532925199433],
#> ID["EPSG",8801]],
#> PARAMETER["Longitude of natural origin",-2,
#> ANGLEUNIT["degree",0.0174532925199433],
#> ID["EPSG",8802]],
#> PARAMETER["Scale factor at natural origin",0.9996012717,
#> SCALEUNIT["unity",1],
#> ID["EPSG",8805]],
#> PARAMETER["False easting",400000,
#> LENGTHUNIT["metre",1],
#> ID["EPSG",8806]],
#> PARAMETER["False northing",-100000,
#> LENGTHUNIT["metre",1],
#> ID["EPSG",8807]]],
#> CS[Cartesian,2],
#> AXIS["(E)",east,
#> ORDER[1],
#> LENGTHUNIT["metre",1]],
#> AXIS["(N)",north,
#> ORDER[2],
#> LENGTHUNIT["metre",1]],
#> USAGE[
#> SCOPE["Engineering survey, topographic mapping."],
#> AREA["United Kingdom (UK) - offshore to boundary of UKCS within 49°45'N to 61°N and 9°W to 2°E; onshore Great Britain (England, Wales and Scotland). Isle of Man onshore."],
#> BBOX[49.75,-9.01,61.01,2.01]],
#> ID["EPSG",27700]]Chapter 8 : Geographic data I and O(Input and Output)
8.1 Retrieving open data
download.file(url = "https://irma.nps.gov/DataStore/DownloadFile/666527",
destfile = "nps_boundary.zip")
unzip(zipfile = "nps_boundary.zip")
usa_parks = read_sf(dsn = "nps_boundary.shp")- 해외여서 접속이 막혀있음
- 공공데이터포털에서 shape 파일 다운받아 불러오기
- 공공데이터포털에서 데이터를 작업 공간에 다운 받기
# unzip(zipfile="C:/202201/GIS/data/부산광역시_교통정보서비스센터 보유 ITS CCTV 현황(SHP)_20210601.zip")
busan <- read_sf(dsn = "./Spatial_Information_Analysis/tl_tracffic_cctv_info.shp", options = "ENCODING:CP949")
busan
#> Simple feature collection with 199 features and 5 fields
#> Geometry type: POINT
#> Dimension: XY
#> Bounding box: xmin: 128.8334 ymin: 35.04753 xmax: 129.24 ymax: 35.35762
#> Geodetic CRS: WGS 84
#> # A tibble: 199 × 6
#> name id url lng lat geometry
#> <chr> <dbl> <chr> <dbl> <dbl> <POINT [°]>
#> 1 센텀중 79 https://its.… 129. 35.2 (129.124 35.17363)
#> 2 수영1호교사거리 53 https://its.… 129. 35.2 (129.1256 35.16539)
#> 3 수영2호교입구(민락) 54 https://its.… 129. 35.2 (129.1281 35.15991)
#> 4 하마정교차로 47 https://its.… 129. 35.2 (129.0678 35.17477)
#> 5 송공(광장) 14 https://its.… 129. 35.2 (129.0688 35.16923)
#> 6 삼전교차로 13 https://its.… 129. 35.2 (129.0638 35.16385)
#> 7 연지삼거리 80 https://its.… 129. 35.2 (129.0553 35.1708)
#> 8 부암교차로 81 https://its.… 129. 35.2 (129.0504 35.16848)
#> 9 동천삼거리 46 https://its.… 129. 35.1 (129.0667 35.13002)
#> 10 진양사거리 82 https://its.… 129. 35.2 (129.05 35.16164)
#> # ℹ 189 more rows
plot(busan)
# unzip(zipfile = "C:/202201/GIS/data/CTPRVN_20220324.zip")
sido <- read_sf(dsn = "./Spatial_Information_Analysis/ctp_rvn.shp", options = "ENCODING:CP949")
sido
#> Simple feature collection with 17 features and 3 fields
#> Geometry type: MULTIPOLYGON
#> Dimension: XY
#> Bounding box: xmin: 746110.3 ymin: 1458754 xmax: 1387950 ymax: 2068444
#> Projected CRS: PCS_ITRF2000_TM
#> # A tibble: 17 × 4
#> CTPRVN_CD CTP_ENG_NM CTP_KOR_NM geometry
#> <chr> <chr> <chr> <MULTIPOLYGON [m]>
#> 1 11 Seoul 서울특별시 (((966987.2 1941111, 966987.1 194…
#> 2 26 Busan 부산광역시 (((1146778 1679624, 1146922 16794…
#> 3 27 Daegu 대구광역시 (((1087860 1760097, 1087860 17600…
#> 4 28 Incheon 인천광역시 (((897904 1961082, 897886.1 19610…
#> 5 29 Gwangju 광주광역시 (((932712.7 1696169, 932781.7 169…
#> 6 30 Daejeon 대전광역시 (((990946.7 1832389, 991057.7 183…
#> 7 31 Ulsan 울산광역시 (((1167950 1710285, 1167950 17102…
#> 8 36 Sejong-si 세종특별자치시 (((971235.9 1844387, 971234.1 184…
#> 9 41 Gyeonggi-do 경기도 (((931607.5 1894480, 931653.8 189…
#> 10 42 Gangwon-do 강원도 (((1163759 1909653, 1163760 19096…
#> 11 43 Chungcheongbuk-do 충청북도 (((1042689 1917663, 1042731 19176…
#> 12 44 Chungcheongnam-do 충청남도 (((862386.1 1805283, 862296.9 180…
#> 13 45 Jeollabuk-do 전라북도 (((902676.5 1717492, 902679.8 171…
#> 14 46 Jellanam-do 전라남도 (((946598.4 1555246, 946595.7 155…
#> 15 47 Gyeongsangbuk-do 경상북도 (((1179681 1750939, 1179685 17509…
#> 16 48 Gyeongsangnam-do 경상남도 (((1053643 1612344, 1053663 16122…
#> 17 50 Jeju-do 제주특별자치도 (((885004.6 1458756, 884996.7 145…
plot(sido)
8.2 지리 데이터 패키지
- rnaturalearth 패키지의
ne_countries()기능을 사용하면 국가 경계 기능을 사용할 수 있음 - osmdata 패키지는 속도가 제한되어 있다는 단점이 있음
- 이러한 한계를 극복하기 위해 osmextract 패키지가 개발
library(rnaturalearth)
#> Support for Spatial objects (`sp`) will be deprecated in {rnaturalearth} and will be removed in a future release of the package. Please use `sf` objects with {rnaturalearth}. For example: `ne_download(returnclass = 'sf')`
usa <- ne_countries(country = "United States of America") # United States borders
#> Warning: The `returnclass` argument of `ne_download()` sp as of rnaturalearth 1.0.0.
#> ℹ Please use `sf` objects with {rnaturalearth}, support for Spatial objects
#> (sp) will be removed in a future release of the package.
class(usa)
#> [1] "SpatialPolygonsDataFrame"
#> attr(,"package")
#> [1] "sp"
usa_sf <- st_as_sf(usa)
plot(usa_sf[1])
korea <- ne_countries(country = "South Korea") # United States borders
class(korea)
#> [1] "SpatialPolygonsDataFrame"
#> attr(,"package")
#> [1] "sp"
korea_sf <- st_as_sf(korea)
plot(korea_sf[1])
8.3 File Formats
- https://r.geocompx.org/read-write.html#file-formats
8.4 Data Input(I)
- gpkg 형식 불러오기
f <- system.file("shapes/world.gpkg", package = "spData")
world = read_sf(f, quiet = TRUE)
tanzania = read_sf(f, query = 'SELECT * FROM world WHERE name_long = "Tanzania"')
tanzania_buf = st_buffer(tanzania, 50000)
tanzania_buf_geom = st_geometry(tanzania_buf)
tanzania_buf_wkt = st_as_text(tanzania_buf_geom)
tanzania_neigh = read_sf(f, wkt_filter = tanzania_buf_wkt)- csv 형식 불러오기
cycle_hire_txt = system.file("misc/cycle_hire_xy.csv", package = "spData")
cycle_hire_xy = read_sf(cycle_hire_txt,
options = c("X_POSSIBLE_NAMES=X", "Y_POSSIBLE_NAMES=Y"))- Well-known text(WKT), Well-known binary(WKB), and the GeoJSON formats
world_txt = system.file("misc/world_wkt.csv", package = "spData")
world_wkt = read_sf(world_txt, options = "GEOM_POSSIBLE_NAMES=WKT")
# the same as
world_wkt2 = st_read(world_txt, options = "GEOM_POSSIBLE_NAMES=WKT",
quiet = TRUE, stringsAsFactors = FALSE, as_tibble = TRUE)- KML file stores geographic information in XML format
u = "https://developers.google.com/kml/documentation/KML_Samples.kml"
download.file(u, "./Spatial_Information_Analysis/KML_Samples.kml")
st_layers("./Spatial_Information_Analysis/KML_Samples.kml")
#> Driver: KML
#> Available layers:
#> layer_name geometry_type features fields crs_name
#> 1 Placemarks 3D Point 3 2 WGS 84
#> 2 Highlighted Icon 3D Point 1 2 WGS 84
#> 3 Paths 3D Line String 6 2 WGS 84
#> 4 Google Campus 3D Polygon 4 2 WGS 84
#> 5 Extruded Polygon 3D Polygon 1 2 WGS 84
#> 6 Absolute and Relative 3D Polygon 4 2 WGS 84
kml = read_sf("./Spatial_Information_Analysis/KML_Samples.kml", layer = "Placemarks")Chapter 9 : Making maps with R
9.1 Static maps
- 정적인 지도는 지리 계산의 가장 일반적인 시각적 출력 유형
plot()또는tmap_mode(plot)
9.1.1 tmap basics
# Add fill layer to nz shape
tm_shape(nz) +
tm_fill()
# Add border layer to nz shape
tm_shape(nz) +
tm_borders()
# Add fill and border layers to nz shape
tm_shape(nz) +
tm_fill() +
tm_borders()
9.1.2 tmap objects
map_nz <- tm_shape(nz) + tm_polygons()
class(map_nz)
#> [1] "tmap"
map_nz
map_nz1 <- map_nz +
tm_shape(nz_elev) + tm_raster(alpha = 0.7)
nz_water <- st_union(nz) %>% st_buffer(22200) %>%
st_cast(to = "LINESTRING")
map_nz2 <- map_nz1 +
tm_shape(nz_water) + tm_lines()
map_nz3 <- map_nz2 +
tm_shape(nz_height) + tm_dots()
tmap_arrange(map_nz1, map_nz2, map_nz3)
alpha: 레이어를 반투명하게 만들기 위해 설정
9.1.3 Aesthetics
ma1 <- tm_shape(nz) + tm_fill(col = "red")
ma2 <- tm_shape(nz) + tm_fill(col = "red", alpha = 0.3)
ma3 <- tm_shape(nz) + tm_borders(col = "blue")
ma4 <- tm_shape(nz) + tm_borders(lwd = 3)
ma5 <- tm_shape(nz) + tm_borders(lty = 2)
ma6 <- tm_shape(nz) + tm_fill(col = "red", alpha = 0.3) +
tm_borders(col = "blue", lwd = 3, lty = 2)
tmap_arrange(ma1, ma2, ma3, ma4, ma5, ma6)
tm_fill()과tm_bubbles()에서 레이어는 기본적으로 회색으로 채워지고tm_lines()은 검은선으로 그려짐- tmap의 인수는 숫자 벡터를 허용하지 않음
plot(st_geometry(nz), col = nz$Land_area) # works
tm_shape(nz) + tm_fill(col = nz$Land_area) # fails
#> Error: Fill argument neither colors nor valid variable name(s)
tm_shape(nz) + tm_fill(col = "Land_area")
- 범례의 제목 설정
legend_title <- expression("Area (km"^2*")")
map_nza <- tm_shape(nz) +
tm_fill(col = "Land_area", title = legend_title) + tm_borders()
map_nza
9.1.4 Color settings
breaks: 색상의 표현 값 범위를 수동으로 설정n: 숫자 변수가 범주화되는 Bin의 수 설정palette: 색 구성표를 정의 (ex.BuGn)
tm1 <- tm_shape(nz) + tm_polygons(col = "Median_income")
breaks = c(0, 3, 4, 5) * 10000
tm2 <- tm_shape(nz) + tm_polygons(col = "Median_income", breaks = breaks)
tm3 <- tm_shape(nz) + tm_polygons(col = "Median_income", n = 10)
tm4 <- tm_shape(nz) + tm_polygons(col = "Median_income", palette = "BuGn")
tmap_arrange(tm1, tm2, tm3, tm4)
tm_shape(nz) + tm_polygons(col = "Median_income", style = "pretty")
tm_shape(nz) + tm_polygons(col = "Median_income", style = "equal")
tm_shape(nz) + tm_polygons(col = "Median_income", style = "quantile")
tm_shape(nz) + tm_polygons(col = "Median_income", style = "jenks")
tm_shape(nz) + tm_polygons(col = "Median_income", style = "cont")
tm_shape(nz) + tm_polygons(col = "Median_income", style = "cat")
style = "pretty": 기본 설정은 가능한 경우 정수로 반올림하고 간격을 균등하게 유지style = "equal": 입력 값을 동일한 범위의 빈으로 나누고 균일한 분포의 변수에 적합(결과 맵이 색상 다양성이 거의 없을 수 있으므로 분포가 치우친 변수에는 권장하지 않음)style = "quantile": 동일한 수의 관찰이 각 범주에 포함되도록 함(빈 범위가 크게 다를 수 있다는 잠재적인 단점이 있음).style = "jenks": 데이터에서 유사한 값의 그룹을 식별하고 범주 간의 차이를 최대화style = "cont": 연속 색상 필드에 많은 색상을 표시하고 연속 래스터에 특히 적합style = "cat": 범주 값을 나타내도록 설계되었으며 각 범주가 고유한 색상을 받도록 함
tm_p1 <- tm_shape(nz) + tm_polygons("Population", palette = "Blues")
tm_p2 <- tm_shape(nz) + tm_polygons("Population", palette = "YlOrBr")
tmap_arrange(tm_p1, tm_p2)
- 순차 팔레트는 단일(ex.
Blues: 밝은 파란색에서 진한 파란색으로 이동) 또는 다중 색상/색조(ex.YlOrBr: 주황색을 통해 밝은 노란색에서 갈색으로 그라데이션)
9.1.5 Layouts
map_nz +
tm_compass(type = "8star", position = c("left", "top")) +
tm_scale_bar(breaks = c(0, 100, 200), text.size = 1)
tm_l1 <- map_nz + tm_layout(title = "New Zealand")
tm_l2 <- map_nz + tm_layout(scale = 5)
tm_l3 <- map_nz + tm_layout(bg.color = "lightblue")
tm_l4 <- map_nz + tm_layout(frame = FALSE)
tmap_arrange(tm_l1, tm_l2, tm_l3, tm_l4)
tm_layout()의 다양한 옵션frame.lwd: 프레임 너비frame.double.line: 이중선 허용 옵션outer.margin,inner.margin: 여백 설정fontface,fontfamily: 글꼴 설정legend.show: 범례 표시 여부legend.position: 범례 위치 변경
tm_s1 <- map_nza + tm_style("bw")
tm_s2 <- map_nza + tm_style("classic")
tm_s3 <- map_nza + tm_style("cobalt")
tm_s4 <- map_nza + tm_style("col_blind")
tmap_arrange(tm_s1, tm_s2, tm_s3, tm_s4)
9.1.6 Faceted maps
urb_1970_2030 <- urban_agglomerations %>%
filter(year %in% c(1970, 1990, 2010, 2030))
tm_shape(world) +
tm_polygons() +
tm_shape(urb_1970_2030) +
tm_symbols(col = "black", border.col = "white", size = "population_millions") +
tm_facets(by = "year", nrow = 2, free.coords = TRUE)
#free.coords : 지도에 자체 경계 상자가 있는지 여부를 지정9.1.7 Inset maps
nz_region <- st_bbox(c(xmin = 1340000, xmax = 1450000,
ymin = 5130000, ymax = 5210000),
crs = st_crs(nz_height)) %>% st_as_sfc()
nz_height_map <- tm_shape(nz_elev, bbox = nz_region) +
tm_raster(style = "cont", palette = "YlGn", legend.show = TRUE) +
tm_shape(nz_height) + tm_symbols(shape = 2, col = "red", size = 1) +
tm_scale_bar(position = c("left", "bottom"))
nz_map <- tm_shape(nz) + tm_polygons() +
tm_shape(nz_height) + tm_symbols(shape = 2, col = "red", size = 0.1) +
tm_shape(nz_region) + tm_borders(lwd = 3)
library(grid)
nz_height_map
print(nz_map, vp = viewport(0.8, 0.27, width = 0.5, height = 0.5))
viewport(): 두개의 맵을 결합
9.2 Animated maps
urb_anim <- tm_shape(world) + tm_polygons() +
tm_shape(urban_agglomerations) + tm_dots(size = "population_millions") +
tm_facets(along = "year", free.coords = FALSE)
tmap_animation(urb_anim, filename = "./Spatial_Information_Analysis/urb_anim.gif", delay = 25)
#> Creating frames
#> =========
#> ====
#> =====
#> ====
#> =====
#> ====
#> =====
#> ====
#> ====
#> =====
#> ====
#> =====
#> ====
#> =====
#> ====
#> =====
#> ====
#>
#> Creating animation
#> Animation saved to D:\Study-Blog\Spatial_Information_Analysis\urb_anim.gifby = year대신along = year을 사용free.coords = FALSE: 각 맵 반복에 대한 맵 범위 유지tmap_animation()을 사용하여.gif로 저장
9.3 Interactive maps
대화형 지도는 데이터 세트를 새로운 차원으로 끌어올릴 수 있음
지도를 기울이고 회전하는 기능과 사용자가 이동 및 확대/축소 할 때 자동으로 업데이트
tmap, mapview, mapdeck, leaflet으로 표현 가능
9.3.1 tmap
tmap_mode("view") #interactive mode
#> tmap mode set to interactive viewing
map_nz
map_nz + tm_basemap(server = "OpenTopoMap")
world_coffee = left_join(world, coffee_data, by = "name_long")
facets = c("coffee_production_2016", "coffee_production_2017")
tm_shape(world_coffee) + tm_polygons(facets) +
tm_facets(nrow = 1, sync = TRUE)tm_basemap()또는tm_options()로 basemap 지정 가능tm_facets()에서sync옵션을 TRUE로 선택하면 여러개의 맵을 동시에 확대/축소할 수 있음
9.3.2 mapview
mapview::mapview(nz)
trails %>%
st_transform(st_crs(franconia)) %>%
st_intersection(franconia[franconia$district == "Oberfranken", ][1]) %>%
st_collection_extract("LINE") %>%
mapview(color = "red", lwd = 3, layer.name = "trails") +
mapview(franconiWa, zcol = "district", burst = TRUE) +
breweries9.3.3 mapdeck
set_token(Sys.getenv("pk.eyJ1IjoiancwMTEyIiwiYSI6ImNsM2ppbzYzNzBrbjQzZHBjMmlocnY2dDUifQ.58-gXpPtvcCGmMt2xEW-ig"))
crash_data = read.csv("https://git.io/geocompr-mapdeck")
crash_data = na.omit(crash_data)
ms = mapdeck_style("dark")
mapdeck(style = ms, pitch = 45, location = c(0, 52), zoom = 4) %>%
add_grid(data = crash_data, lat = "lat", lon = "lng", cell_size = 1000,
elevation_scale = 50, layer_id = "grid_layer",
colour_range = viridisLite::plasma(6))
#> Registered S3 method overwritten by 'jsonify':
#> method from
#> print.json jsonlite